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ABSTRACT 



We have developed a three dimensional Monte Carlo photoionization code tailored 
for the study of Galactic H ii regions and the percolation of ionizing photons in 
diffuse ionized gas. We describe the code, our calculation of photoionization, heating 
& cooling, and the approximations we have employed for the low density H ii regions 
we wish to study. Our code gives results in agreement with the Lexington H ii region 
benchmarks. We show an example of a 2D shadowed region and point out the very 
significant effect that diffuse radiation produced by recombinations of helium has on 
the temperature within the shadow. 
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1 INTRODUCTION 

A major characteristic of the interstellar medium (ISM) 
is that it is clumpy, with filamentary structure over 
scales ranging from hundreds of parsecs through sub-AU 
cloudlets. The spatial density of electrons is approximately 
a power law over an a stonishing 12 orders of magnitude 
iArmstrong et al lEiii), from ~300 pc to <0.1 AU in the 
ISM. The magnetic field of the Galaxy varies in a simi- 
lar fashion on scales from 0.01 to 100 pc. Molecular clouds 
show a fr actal (scale-free) spatial variation of the intensity of 
CO fe.g..lFalgarone. Ph illips, fc WalkeJlQQlT) and neutral H 
i|st animirovic fc Lazarian>.200Ll . Sub-parsec variations are 
shown by the variation of ionic column densities along dif- 
ferent ISM sightlines towar ds individual stars in close bi- 
naries iLauroesch. Mever. fc B lades 2000J. L arge area ve- 
locity resolved surve ys of the molecular (e.g.. [ Simon et all 
l200l]f) neutral (e.g.. iHartmann fc Burtqn|[l99ar and ion- 
ized iRevnolds. Haffner^^Mad^er3l200^) gas now allow us 
to probe the geometry, kinematics, ionization, and temper- 
ature structure of the ISM. 

Our development of Monte Carlo radiation transfer 
techniques to study photon scattering, ionization, and ra- 
diative equilibrium dust temperatures in complex geome- 
tries is motivated by these and the many other observations 
that clearly require three dimensional modeling analyses. 
In this paper we describe improveme nts and extensions to 
our Monte Carlo photoionization code JWood fc Loebl2000l) 



that will enable us to model the 3D temperature, ioniza- 
tion structure, emission line images, and spectra of low den- 
sity H II regions and the Diffuse Ionized Gas (DIG), often 
called "the warm ionized medium" or "Reynolds layer" . This 
diffuse ionized gas has a complicated structure and surely 
requires three-dimensional modeling, both as regards the 
sources of its ionization and the transfer of the ionizing ra- 
diation. How can Lyman continuum photons from O and B 
stars percolate from their midplane origins to ionize the high 
latitude gas? A non-uniform ISM is required that provides 
low density paths for ionizing photons to reach the halo. Us- 
ing a "Stromgren volume" io nization technique an d a sta- 
tistical approach to clumping, iMiller fc Cox l flQOS") showed 
that O and B stars can ioniz e the D IG (see also Dove & ShulJ 
ir994l:lDove. Shull. fc Ferraral2000l). Using Mon t e Carl o sim- 
ulations for hvdrogen ionization. ICiardi et alJ (|20o3) have 
investiga ted ionization str uctures in a fractal ISM as pro- 
posed bv lElmegreerJ 1^93)- In Elmegreen's picture, the ob- 
served Ha emission may arise from the ionized surfaces of 
fractal clouds. This suggestion clearly needs to be investi- 
gated with 3D radiation transfer techniques. 

Other features of the DIG that may require 3D 
analyses are the temperature structure and relative ion- 
ization structure of He and H. The temperature of 
the DIG, as determined from l ine intensity ratios (e.g., 
iRevnolds. Haffner. fc TuftdllQQgTl . appears to require addi- 
tional heating over and above that provided by photoion- 
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ization heating llMathij|200(fl . Among the possible heating 
mechanisms, dissip ation of energy in a tur bulent medium 
has been proposed (|Minter fc SDa ngler"l99(fl. definitely 3D! 
Another question is, can clumping /shad owing explain the 
appa rent lack of He" in H'*' regions jReynolds fc Tuftd 
Il99fj) ? In this scenario, direct stellar photons are prevented 
from ionizing He in shadow regions behind dense clumps, 
but diffuse ionizing photons from H and He recombinations 
can percolate into the shadow regions and ionize H, but not 
He. 

Our work on the analysis of scattere d light in fractal 
reflection nebulae ijMathis, Wh itney, fc W ood 2002) high- 
lighted the crucial role of geometry when trying to de- 
termine dust properties from scattered light observations. 
Much of what we know about abundances of elements 
comes from analysis of spectra utilizing one-dimensional 
photoionization techniques (see e.g. ^Stasihska 2002). Over 
the last decade several three dimensional photoioniza- 
tion codes have been developed. T he code presented by 
iBaesseen. Diesch. fc Grewind (Il990l) uses the "on the spot" 
approximation for the diffuse radiation field, while that of 
[Gruenwald, Viegas. & BroEuicrc ( 1993) employs an adap- 
tive mesh grid and an iterative procedure to determine the 
contribution of the diffuse radiation field. Monte Carlo tech- 
niques naturally include diffuse radiation within three di- 
mensional systems and Monte Carlo photoionization codes 
have been described by .Ercolano et al. ( 2 003') for the three- 
dimensional case and lOch. Lucv. fc Rosalll998il for the one- 
dimensional case. 

We have extende d our basic hydrogen- only Monte Carlo 
photoionization code iWood fc Loebl2 000') to treat the pho- 
toionization, heating, & cooling for the complex geometry 
expected in low density H ll regions a nd the diffuse ISM. 
Our co de i s similar to those descr ibed bv lOch. Lucv. fc Rosal 
Jl998h and lErcolano et alJ (|200!i), but calculates the rates of 
photoionization and heating "on the fly" , rather than calcu- 
lating the mean intensity of the radiation field (see also Lucy 
1999). It also treats reproce ssing of radiation in a diff erent 
manner to that described bv lOch. Lucv. fc Rosal lll998ll . us- 
ing "photon packets" rather than "energy packets." 

In §2 we describe our assumptions and approximations 
to model the DIG. §3 describes the Monte Carlo technique, 
our method of discretizing the radiation field into "photon 
packets" rather than "energy packets," and how we deter- 
mine the diffuse radiation field using ratios of recombination 
rates. §4, §5, fc §6 deal with our photoionization and heat- 
ing/cooling algorithms, §7 presents the atomic data used in 
our code, and §8 compares our code with standard H II re- 
gion benchmarks and gives an example of a 2D calculation. 



2 APPROXIMATIONS FOR THE DIG 

We make the following assumptions and approximations for 
the opacity, heating and cooling. 

• We consider only radiation energetic enough to ionize 
H (hu ^13.6 eV). Photons of lower energies are ignored un- 
til we calculate the emergent intensities of certain spectral 
features after the simulation is completed. 

• All ions are in the ground state when they absorb ra- 
diation or recombine with an electron. This is the standard 
"nebular approximation." The low-lying levels giving rise to 



fine structure lines might be an exception. Radiative transfer 
in these lines would require consideration of the populations 
of the levels involved in their production. We assume that 
the fine structure lines are optically thin. 

• No H-ionizing photons are produced by recombinations 
to excited levels of H or He. This is because fcTe ^ x(H'') 
= 13.6 eV, the ionization potential of H. 

• There is no coUisional ionization or de-excitation. In 
effect, this limits the densities we can study to be less than a 
few X 10^ cm~'^, which is generally above that encountered in 
the Diffuse Ionized Gas and typical H II regions, but within 
the typical range of planetary nebulae. 

• There is no He^^ because the stellar radiation field has 
negligible luminosity above 54.4 eV, consistent with obser- 
vations of almost all H ll regions. 

• The opacity is due to continuous opacity from neutral H 
and He; we ignore contributions from other ions because of 
their low abundances. At present we ignore all other opacity 
sources (e.g., dust, UV metal line blanketing, etc.) except 
photoionization from oxygen. 

• Photoionization and heating by the He I Lyman a line 
is included in an "on-the-spot" approximation (see >l4.1l for 
a discussion of the physics). 

• We keep track of the abundances of the stages of ioniza- 
tion of H, He, C, N, O, Ne, and S that have their ionization 
potentials ^ x(H°) and < x(He+) = 54.4 eV. These are H°, 
He°, C+ - C+^ N° - N+^ 0° - 0+^ Ne° - Ne+^ and S+ 
- S"""^. These chemical elements are chosen because some of 
their ions are important for the heating and cooling rates. 
It is straightforward to include other elements, but they do 
not appreciably affect the H and He ionization structure or 
contribute to the heating or cooling. We assume that there 
is no S" because its ionization potential is 10.36 eV, so it is 
ionized by the ambient stellar radiation fields softer than hv 
= 13.6 eV (see Savage & Scmbach 1996 for observations). 

• Heating arises solely from the photoionization of H and 
He. Cooling is from recombinations of H and He, free-free 
emission, and coUisionally excited line radiation from N, O, 
Ne, S (see below). 

• We ignore dynamics, turbulence, and the effects of 
shocks and ionization fronts. 

It is straightforward to include more ions and the effects 
of dust opacity. The assumption of no opacity from heavy 
elements may also be easily relaxed. 



3 MONTE CARLO PHOTOIONIZATION 

We want to determine the time-independent ionization and 
temperature structure of a nebula with an assumed fixed 
density distribution. We do this by Monte Carlo radiative 
transfer followed by iterations to find the ionization struc- 
ture that determines the opacity. 

Within each iteration we consider packets of ionizing 
photons emitted from the exciting star(s) of the nebula. We 
choose the frequency from the probability distribution of 
the photon emission from the sources. We follow each stellar 
photon packet through the density grid until it is absorbed 
and causes photoionization. The subsequent recombination 
often results in the production of another ionizing photon 
packet, thereby forming the diffuse radiation field. The fre- 
quency of the new packet is sampled from various atomic 
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probability distributions. We first determine if tfie packet 
was absorbed by H° or He", and then determine tlie fre- 
quency of tfie re-emitted packet (see below) . We continue the 
Monte Carlo random walk of absorption/re-emission events 
until the packet is re-emitted with hi/ < 13.6 eV, at which 
time we terminate it. 



3.1 Photon Packets versus Energy Packets 

Our approach diff ers somewhat from pre v ious Monte 
Carlo simulations l Och. Lucv. fc Rosal Il998t ILucvI Il999l: 
iBiorkman fcW ood 200ll: lErcolano et alj2o'oii that followed 
the energy of a packet, rather than its number of photons. 
Considering energy is natural for investigating problems of 
radiative equilibrium (e.g., scattering and thermal emission 
by dust, or stellar atmospheres), since the energy that is ab- 
sorbed is thermally re-emitted with exactly the same energy 
but with a source function, S^, given by Kirchhoff's Law 
{Si, = jv/kv = Bu(T), where ji, is the emissivity and kv the 
absorption coefficient). Kirchhoff's Law is applicable only 
when the populations of atomic levels are in Local Thermo- 
dynamic Equilibrium and does not apply to nebulae. The 
number of photons associated with an energy packet varies 
with its frequency, so it is not as straightforward to apply 
atomic probabilities in following the fate of an energy packet. 
The energy of an energy packet is given by e = L At/N and 
is kept constant for radiative equilibrium codes, where L is 
the stellar luminosity. At is the time covered by the simula- 
tion, and A'^ is the number of packets used in the simulation. 
The energies of our packets are e = QhvAt/N, where Q 
is the ionizing photon luminosity (units s~^). The energy 
is not kept constant as the packet changes frequency, but 
keeps its number of photons constant (see also Lucy 2003). 

Since our initial goal is to determine the 3D ionization 
and temperature structure, we do not follow the non-ionizing 
diffuse photons: they are emitted with hv < 13.6 eV and will 
escape from the nebula, so we terminate the ionizing photon 
and start a new one from the star. We thus save memory 
by not having to form probability distribution functions in 
each cell for finding the frequency of non-ionizing photons. 
Rather, we determine a converged ionization and tempera- 
ture structure and then employ a second code to integrate 
numerically through our resulting emissivity grid to form 
the emergent spectrum and images in various spectral lines. 

Current desktop machines using 1Gb RAM can handle 
grids of around 128^ cells for our code. The iterative pho- 
toionization and heating calculation (described in §4 and 5) 
updates the temperature and ionic fractions of all species 
we are tracking in each cell throughout the grid. At the end 
of each iteration we have a new ionization structure and 
opacity grid for the next iteration. 



3.2 Emission Direction and Frequency 

In Monte Carlo techniques we must decide on the outcome of 
events that depend on some physical variable, say x (e.g., the 
frequency v), with a given probability, P{x). The important 
quantity is the cumulative distribution function, C{x), given 

by 



C{x)^ 



r.n^Qdx' 

/;+ P(rc')dx' 



where x- is the minimum value of x and x+ the maximum. 
We find a randomly selected value of x by casting a random 
number'^, ^, in the range [0,1] and inverting the relation 
^ = C{x). 

We assume the density structure of our nebular simula- 
tion, possibly quite complex as appropriate to the ISM. We 
begin the Monte Carlo radiation transfer by emitting pho- 
ton packets from the exciting star(s) with an initial guess for 
the ionization structure. The initial guess for the ionization 
structure does not affect the final converged solution, but 
it does affect the number of iterations required for conver- 
gence. For isotropic emission, the random directions are 



cos 6* = 2^ - 1 
= 2^^ , 



(1) 
(2) 



where 9 and are the usual spherical polar coordinates. The 
direction cosines for the initial photon trajectories are then 
given by 



Tlx — sin 9 cos ( 
Uy — sin 9 sin q 
n, ~ cos S ■ 



(3) 



The photon frequency must be sampled to reproduce 
the ionizing spectrum of the source, L{v), which may be a 
blackbody or model atmosphere. We pre-tabulate the cumu- 
lative photon distribution function. 



c(i^O 



r' L(v)/huAv 
r°° L(v)/hvdu 



(4) 



A random frequency is chosen by solving for v in ^ — C{u) 
by interpolating in the pre-tabulated cumulative distribu- 
tion function C{vi). From v we calculate the photoionization 
cross sections av{E^) and 0^(116") from analytic formulae. 
The optical depth, r^, along a path length 5* is 



[n(H°)a4H°) + n(He'^)a,(He'')] ds 



(5) 



where n(H°) and n(He°) vary from cell to cell. We generate a 
random optical depth by r = — log ^ and integrate through 
the opacity grid until the distance 5* has been reached. If 
the distance S lies outside the simulation grid the photon 
has exited our simulation and we start another packet from 
the star. 



3.3 The Diffuse Radiation Field 

On finding the position associated with the randomly cho- 
sen optical depth, we decide if the photon packet is ab- 
sorbed by H" or He'^. The packet is absorbed by H" if 
^ ^ n(H°)a^(H°)/[n(H'')a^(H°) -f n(He°)a^(He°)], and by 

^ A new r andom number is generated by the subroutine 
ran2(iseed) jPress et abligO^) each time ^ appears in this paper. 
A new integer, iseed, is automatically generated by the subroutine 
after each entry. 
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He'' otherwise. If absorbed by H", it will be re-emitted either 
as a Lyman continuum photon (and subsequently tracked 
with the Monte Carlo technique) or terminated. If absorbed 
by He" there are several energy channels into which it can 
be re-emitted as described below. Once a photon packet has 
been absorbed we re-emit it with a frequency chosen using 
the following algorithm. 

3.3.1 HI Emission 

The probability that the photon is emitted in the hydrogen 
Lyman continuum is 



P(H I Ly c) = 



Qi(H°,re) 
OA (HO ,T, 



(6) 



where ai is the recombination coefficient to the ground level 
and a A the coefficient to all levels, including n = 1. The 
packet emerges in the Lyman continuum if f f (Ly c). If 
the packet is chosen to be re-emitted in the Lyman contin- 
uum, we sample its frequency using the C{u) derived from 
the Saha-Milne emissivity: 



CLyc( 



. 27r me fc T, 

X a.(H")e-'^(''-''o)/fcTe 



n 



;H+)ne 



(7) 



(8) 



We pre-tabulate Ci^y^iv, Te) and from ^ = C(y,Te) we 
determine v from ^ by interpolating in the table. The packet 
is now re-emitted with a direction chosen from an isotropic 
distribution (Eqs. 1 and 2). 

At the start of each iteration we form a 3D grid that 
contains P(H I Ly c) and use this grid to decide the fate of 
photon packets absorbed by H''. Similar grids are formed for 
the He emission probabilities as described in the next sec- 
tion. These probability grids are updated after each iteration 
since the recombination rates are temperature dependent. 



3.3.2 He I Emission 

An absorption by He is followed by recombination and pos- 
sibly by radiative cascades into either the metastable levels 
2^S or 2^S, to 2^P (the upper level of Lyman q), or to the 
ground level l^S. The probabilities for emission in these four 
channels are as follows: 

(i) He I Lyman continuum {hf > 24.6 eV) is analogous 
to that given for H in 

(ii) The probability of emission as a 19.8 eV photon 
packet from the transition 2^S ^ l^S is 



P(19.8eV) = 



a°fg(He°,r,) 
aA(HeO,Te) 



(9) 



The effective recombination rate takes into account all 
means of populating the given level, both direct recombi- 
nation and via cascades from higher levels, from which low 
energy photons are emitted. 

(iii) The probability of emission from the He I two-photon 
continuum from 2^5* is 



P(He I2q) = 



QA(HeO,r,) 



(10) 



Note that a"^ gCHe'^ ,Te) does not include de-excitations 
from 2^P, which will be discussed in §4.1 

A fraction 0.28 of the photons in this emission channel 
is able to ionize H", so the number of photons emitted per 
decay is 2(0.28) = 0.56 photons. After determining that the 
packet is to be emitted in the 2-photon continuum, we see if 
^ 0.56. If so, we choose a freq uency from a pre-tabulated 
Civ) for two-photon emission jPrake. Victor, fc Dalgarnd 
IT969L Table II) . 

(iv) The probability of emission as a He I Lyman a pho- 
ton (hu = 21.2 eV) is 



P(HelLya) 



Q^2?p(He°,Te) 



(11) 



QA(HeO,Te) 

lOch. Lucv. fc Rosal Jl998l) and lErcolano et al.l J2003h fol- 
lowed energy packets of several of the He I Lyman lines, as- 
suming that the various He Lyman lines are absorbed only 
by the H° continuous opacity, propagating through the grid 
and maintaining their identity without degrading into Ly- 
man Q. We do not follow this prescription because resonance 
lines will be scattered many times without diffusing far in 
space. Following each scattering from n}P {n ^ 2), there is 
a finite probability of re-emission into a lower excited level 
rather than to 1^5*. Subsequent emissions from that excited 
level cascade downwards and populate either 2^5* or 2^ P. 
Conversion to 2^S is followed by He I 2-photon continuum 
emission. We discuss our treatment of the 2^P in 

To form the various probabilities, we use analytic fits to 
the temperature dependent recombination rates for H and 
He. These rates, along with all other atomic data, are pre- 
sented in J7| Since there are only four He I levels with ap- 
preciable populations, we have 



I cff I cff I cff 



(12) 



We also assume that recombinations to triplets produce only 
the 19.8 eV 2 ^S ^l^S t ransition. We do not follow 2='P 
l^S, as did Oc h. Lucv. fc Rosa ( 1998), since at low densities 
2'^P converts to 2^5* via 10830A emission. 



4 PHOTOIONIZATION 

We must calculate the photoionization, heating, and cooling 
in order to update the ionic fractions and electron temper- 
ature after each iteration. 

The basic ph otoionization equations, in the notation of 
iQsterbrockl (ll989D . are 



4nJ^ 
hv 



a^{X+')Au = n{X+'+'-)n,a{X+\T,), (13) 



where n(X^') and n{X^^^^) are the number densities (per 
cm"') of successive stages of ionization, Vi is the threshold for 
photoionization from Ue is the electron number density, 
and a{X^^ ,Te) is the recombination coefficient to all levels 
of X^^ . In our code the photon luminosity of the sources, 
Q, is divided into A'' packets, so we calculate the integral on 
the left of eqn. 11311 by 



4:11 J ^ 

hv 



N AV 



^la,{X^ 



(14) 
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where the summation is over all path lengths through the 
volume /S.V for photons in the frequency range [v, v + Av\ 
and I is the path length of each photon counted in the sum 
(Lucy 1999, 2003). 

Thus, we do not explicitly calculate the mean intensity 
throughout the grid in order to compute the integral of 
^'Kj,,ai,(X^^^lhv for each ion considered. To save J^, we 
would need storage for each grid cell with as many frequen- 
cies as being considered. 

We add the contribution of each ionizing photon packet, 
regardless of stellar or nebular origin, to the sum of 
?ai/(X+') maintained for each ion, within each cell in the 
grid. At the end of the iteration we have a 3D array for each 
ion that contains the integral (x, y, z). The photoioniza- 
tion rate is then n(X+*) l^+i within each cell. So long as the 
number of ions being tracked is less than the typical num- 
ber of frequencies required for an accurate representation of 
the mean intensity, our technique will save significantly on 
memory storage. We use the above equations to solve the 
ionization balance of most ions that we follow. Exceptions 
to this are our treatment of He I Lya photons that can ion- 
ize H", and also charge exchange that couples the ionization 
of some ions to H and He. 

We now outline our treatments of the ionization balance 
for these situations. 



4.1 He I Lyman a 

We do not need to follow the diffusion of the He I Lya 
resonance line {2^P l^S) in space and frequency. The 
level 2^ P is depopulated either by emission of a He I Lya 
photon {2^P l^S) or a 2.06^m photon (2^P 2^S). 
Due to the large resonant line opacity, the He I Lya will be 
scattered locally many times without traveling far. It will 
eventually be absorbed locally ( "on the spot" ) by H" or else 
decay to 2^S. 

Note that the "on the spot" approximation is used due 
to the large He I Lya resonant opacity and not the H'' opac- 
ity. 

What is P(HoTs), the probability of the He I Lya 
photon being absorbed by H" as opposed to the transition 
2^P — > 2^S occurring? The average number of resonant 
scatterings before decay to 2^S is A pip -,iq /Aqi p o = 914 
iBeniamin. SkiUman. fc Smitd Il999t) . The mean free path 
between scatterings, /scat, is [n(He°) a^(Lya)]~^, where the 
line opacity is 



a^(Lya) = {tt^ e'^ /m^c) /21 exp[~{Au/AuD) ]/AuD , (15) 

where /21, the oscillator strength of the transition, is 0.29. 
Af/Ai/D is the displacement of the frequency of the pho- 
ton from the line center, in Doppler widths. This displace- 
ment is determined by the speed the atom had before re- 
combination, and so has a Maxwellian probability distri- 
bution (if we neglect any non-thermal components in the 
local velocity distribution) . The Maxwell-Boltzmann distri- 
bution can be written fMB{v) oc exp[—{v/vth)^] and 
Av/AvD = (w/uth), with «th = {2kT/My^^. We average 
over the line profile: 



cxp[-(tj/Dth)^l cxp[-(u/jJti,)^l 
f°° j)2 cxp[-(tj/Dth)^] dv 

u 

= 2-^/2 . (17) 

The number of absorptions by H° during the 914 scat- 
terings is 914 /scat n(H°) a^(H°), with a^^°) = 1.9 x 10"^** 
cm^. The probability of the photon being absorbed by H'' is 

n(HO)a.(HO) 



P(HoTs) [„(HO)a^(HO) -Hn(He")a,(HeO)/914] 

= (l + 0.77/(He°)/[/(H«)(T/10*K)^''^])~' , (18) 

where /(He") and /(H*^) are the neutral fractions of He and 
H, respectively, and He/H = 0.1 in abundance has been as- 
sumed. The /(He°)//(H°) favors He, especially for softer 
exciting stars. For instance, the inner parts of the nebula 
surrounding a 40000 K black body favor He" by a factor 
of ~7, so 80% of transitions to 2^P are converted to 2}S 
rather than being absorbed by H". For harder exciting spec- 
tra, the /(He")//(H'') is <2, and a significant number of 
Lya photons are absorbed by H". Heating by Lya is signif- 
icantly greater than from 2-photon emission. However, the 
(2^5' ^ l^S) photons are far more important than any from 
the singlets as regards both ionization and heating. 

In the "on-the-spot" treatment, local absorptions 
of Lya are exactly balanced by local emissions. The 
rate of "on-the-spot" absorptions by H" is P(Hots) 
times the rate of recombinations into He°(2^P), or 
P(HoTs)n(He+)nca'''*(He°,2ip). These absorptions by H" 
are also given by ?i(H") ai/(H")[47r J^(Lya)//ii^]. By equating 
these rates, we obtain a modified ionization balance equation 
for H: 



n(H°) /_J(47rJ„//ii/)a^(H°)di^ = 
,[n(H+)a(H",Ve) - P(HoTs)n(He+)a°?p] , 



(19) 



{exp[-(A/./Ai/i5)2]) 



(16) 



where the does not include He I Lya, which is contained 
in the final term on the right side of the equation. 

4.2 Charge Exchange 

The rate of exchange of an electron between an ion (X^*) 
and H° or He" can be comparable to or exceed the radia- 
tive recombina tion rates. With our adopt ed charge exchange 
cross sections llKingdon fc Ferlandlll99B) . the rate of transi- 
tions into N^^, O", and O"*" are seriously affected or dom- 
inated by exchange with H". The opposite effect, in which 
charge exchange dominates the ionization of the lower stage, 
takes place for 0°. A result is that in regions where 0° is 
important (i.e., H" is appreciable), (0°/0+) ~ 9/8 (HVH+) 
for all To because t he ioni zation potential of O'' is almost the 
same as for H'' fsee lOsterbrock 1989 . p. 42 for a discussion.) 
For N", charge exchange is not dominant but not negligible. 
We included the reactions in the rates for reactions resulting 
in and as well. 

With charge exchange included, the (O^/O^) equation 
becomes 

n(0°) K / 47r J,//ii/d!/ + n(H+)Ccxch(0°, H+)] = 
n(0+) K arcc(0+) + 7i(H0) Ccxch(0+, H")] . 
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The situation regarding charge exchange with He'' is less 
clear because the cross section is very difficult to compute. 
We account for charge exchange with He" for interactions 
resulting in N+, N+^, and 0+. 



5 HEATING & COOLING 

The temperature is determined by finding the root of 
G(r) = C(T), where G and C are the heating and cooling 
rates due to the p rocesses described below (see Chapter 3 of 
IOsterbrocklll989li . 



5.1 Heating 

We consider only heating by photoionization of H and He 
and calculate the heating rate in each grid cell in the same 
manner as we calculated the photoionization rate. For ex- 
ample, for photoionization of He the heating rate is 



G(He'') = 



47rJ^ 
hv 



a^(He'') hiu — v^^o) dv , 



(20) 



so we form integrals analogous to the photoionization rate 



[Q/{N'AV)]J2laA^e°){hv - huti.o) . 



(21) 



As with photoionization, we assume "on-the- 
spot" heating by He I Lya packets, so a term 
P{UoTs)n{Re+)n^al^p{Lya)[hu{Lya) ~ x(H'')] is added 
to the heating per volume of each cell. 



5.2 Cooling from Recombination and Free-Free 
Radiation 

Cooling rates as a function of temperature, C, due to recom- 
bination of H^ and He""" are taken from Hummer (1994) and 
Hummer & Storey ( 1998). The cooling rate due to free-free 
radiation is given bv lOsterbrockl 1^83): 

Cff = 1.42 X IQ-^'^ [n(H+) + n(He+)] rie gaVfe , (22) 

where we use the Katz ot "ai] \l99di) fit to the free-free Gaunt 
factor data in Soitzcr ( 1978) 

gs = 1.1 + 0.34 exp [-(5.5 - logr,)V3] . (23) 



5.3 Other Cooling 

We include coUisional excitation of the lowest five levels of 
the ions of N°, N+,0^0+,0+^ Ne+^ S+, and S+^ and of 
2 lev els of N^^ and Ne"*", as in other photoionization codes 
fe.g.. lErcolano et al.ll2003 ). For this cooling by coUisionally 
excited forbidden lines we use the Einstein A values and 
collisio n strengths from the compilation by Pra dhan fc Pend 
(^2^ as updated on A. Pradhan's webpage.^ 



' |http : / / www-astronomy. mps . ohio-state. edu/~pradhan/ 1 



6 SOME CODE DETAILS 

After each iteration, we simultaneously solve for the ioniza- 
tion and temperature using the equations outlined above. 
The code will converge quickly if we have a good first guess 
of the ionization and temperature structure. However, for 
the 3D problems that we wish to address, such a guess will 
be very difficult to make. We therefore start the first itera- 
tion assuming that H and He are fully ionized, so all stellar 
photons exit the grid on the first iteration, but do contribute 
to the mean intensity counters. If we started with a neutral 
grid, photon packets could not diffuse far from their point 
of origin due to the large neutral opacity and convergence 
would be very slow, if at all. 

As in lErcolano et al.l tOQ^ . we run several iterations 
with a small number of photon packets to get a rough idea 
of the ionization and temperature structure: 10^ packets 
for the first eight iterations, up to 10* for higher iterations 
^20. Our models are well converged within 12 iterations; 
the higher iterations with large numbers of photon packets 
provide higher signal-to-noise for the ionization and temper- 
ature structure. 

In 3D simulations, the diffuse field ionizes walls strongly 
inclined to the direction to the star. We find a good early ap- 
proximation to the diffuse field by initially setting all cells to 
be fully ionized, thereby predicting the diffuse radiation field 
if the nebula were fully ionized. For the next iteration, if the 
stellar radiation cannot maintain a high degree of ionization 
the outer material becomes almost neutral, since the diffuse 
field is much weaker than the stellar field. The structure 
reaches a good approximation to the eventual diffuse field 
after about five iterations. If the shadowed regions were al- 
lowed to turn neutral before the diffuse field had converged, 
we would face the same convergence problem (with each new 
iteration ionizing a small zone of previously neutral mate- 
rial) described above. 

We run our code on a single desktop mach ine within a 
reasonable CPU time (see lErcol ano et'al ] |2003l for a discus- 
sion of code speed up via parallelization) , running ~1.0M 
photon packets per minute on a 2.4 GHz processor. 



7 ATOMIC DATA 

We use standard references for the line and continuum 



emissivities (iBrocklehurslJ Il972': 


Brown & Mathews 


'l970|; 


IStorev & HummeJ 


Il995t iBeniamin. Skillman. & SmitsI 



the fitting formula of IVerner. Ferland. Korista. fc YakovlevI 
^96. The recombination rat e to all levels aA(H°,T 'e) is 
given by the fitting formula in lVerner fc Ferlandl il996l) . For 
the recombinations to the ground state, l^S, we fitted the 
temperature dependent rates in Osterbrock (1989, Table 
2.1), with 

ai{}f,T,) = 1.58 X lQ-'\n/lQ'^K)-°-^^ . (24) 

For He we use the lVerner. Ferland. Korista. fc YakovlevI 
l|19963 formula for the photoionization cross section. For the 
re combination rates we use the fitting formulae in Table 1 
of IBeniamin. Skillman. fc SmitsI il999l) that give us the di- 
rect recombination rates to the ground level, l^S, and the 
effective rate to 2^5*. We calculated the effective rates to 
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2^P and to 2^5, ignoring the 2^P 2^S transition be- 
cause we separately consider its competition with absorp- 
tion by H (see the discussion in §4.1). The rates to 2^P 
and 2^S were formed by adding the dire ct recombination 
rates iBeniamin. Skillman. fc Smitj Jl999ft . Table 1), effec- 
tive rates for n > 5 (Table 3), and effective rates from all 
lines for 3 ^ n ^ 5 that connect to the appropriate lower 
level (Table 5). In summary, 

ai(He",r,) = 1.54 x 10~"(re/10*)-°-^*^ 
a^T5(He",re) = 2.06 x 10-"(re/10'*)-°-^^® 
Q2ip(He°,re) = 4.17 X 10-"(re/10*)-''-*" 
a^fs(He^re) = 2.10 x 10""(re/10*)-° ''^'* , (25) 

thus enabling us to form the probabilities for reprocessing 
packets by He ( H3.3.2II . 

The ionization cross section for the ions we use are taken 
from IVerner. Ferland. Korista. fc Yakovlev| l|l^ ^^^ and the 
radiat ive recombination rates are from IVerner fc Ferlandl 
il996h . Also important in the recombination are dielectronic 
recombination rates, for which we use values determined by 
Nussbaumer & Storey (1983, 1984, 1987) for N, O, and Ne. 
The dielectronic recombination rates for sulphur have not 
been calculated, and we normally use those recommended 
bv lAli e t al. (1991) for the first four ionization stages of S: 
3 X 10"", 3 X 10"^^ 1.5 X 10~", and 2.5 x 10"\ which are 
the mean ra tes for the first fo ur atoms/ions of C, N, and O as 
discussed b vlAli et al.ll^l991^. The cha rge exchange rates are 
taken from Kinadon & Forland (199?) for reactions with H 
and from .Arnaud fc R.othenfiue. (1985" ) for those with He". 



8 SOME RESULTS AND BENCHMARKS 

8.1 Lexington H II Region Benchmarks 

Any 3D photoionization code must reproduce two bench- 
mark tests of blackbody exciting stars ionizing uni- 
form spherical nebulae of a specified composition. The 
benchmarks have been discussed extensiv ely in workshops 
JPerland et al.lll995llP6auignot et ai]l200ll: see these for the 
details of the input parameters). Results of many codes are 
given in Table 1 (for a 40kK central star) and Table 2 (T, = 
20kK). The authors of the other codes are given in the foot- 
note to Table 1. The dielectronic recombin ation coefficients 
assumed for S are zero for and those in lAli et al.1 (1199111 
for and We believe these were adopted by most of 
the models in each table. Some models in the benchmark 
took the optical depths of fine structure lines into account, 
but we did not. Our results are given in the last column in 
each table. 

Figure 1 shows plots of ionization and temperature with 
radius for the =40000 K benchmark and compares our 
results (cro sses) with the Monte Carlo photoionization code 
MOCASSIN iErcolanp et alJl20 03!l (squares) and the widely 
used code cloudy i^Handl l994') (lines) . Figure 2 presents 
the same quantities for the T* = 20000 K benchmark. In 
comparing with MOCASSIN, we set up the density grids to 
have the same spatial resolution: our code used a grid of 
65^ cells and mocassin, which utilizes the symmetry of the 
benchmark, performed the radiation transfer in an 33^ oc- 
tant. For these simulations our code used 10* photon pack- 



ets and mocassin used 6.4 x 10^ energy packets in the fi- 
nal (20th) iteration, though both codes had essentially con- 
verged by iteration 14. Among the three codes, the ioniza- 
tion fractions are all in very good agreement. Because of our 
65^ grid resolution, our code produces slightly more ioniza- 
tion at the outer edge. There is a larger spread in the pre- 
dicted temperatures for the hot benchmark, with our code 
being slightly hotter in the inner portions of the nebula and 
MOCASSIN being slightly cooler towards the edge of the neb- 
ula. The agreement for the cool benchmark temperatures is 
very good. 

The large spread of some ionization fractions (e.g., 
0'''^/0 in Fig. 2) illustrates a property of Monte Carlo codes 
unless special steps are taken. The ionization potential of O"*" 
is 35.1 eV; an ionizing photon is so rare that even with 10* 
photon packets there is considerable statistical noise. If an 
accurate prediction of the (O^^/O) were one of our goals, 
we could have forced energetic stellar photons to be emitted 
and then correct the results by the known probabilities of 
the events we forced. 

The cubic nature and limited spatial resolution of our 
grid makes the comparison problematic for one ion, 0°, 
mainly found in the very outermost regions of the nebula. 
In these regions, {O^ /O'^) = 9/8 H^/H"*" because of a reso- 
nant charge exchange between O"^ and H''. There is a rather 
sharp peak in the temperature because of the hardening of 
stellar radiation caused by selective absorption within the 
nebula, and a steep temperature gradient because of the 
large opacity. With better spatial resolution, the volume oc- 
cupied by one cell would contain regions of almost neutral 
hydrogen, so the temperature assigned to the whole cell, 
and the strength of [O I] A6300, are of limited meaning. We 
present results from our cells with H'^/H ^ 0.25; the cells 
more neutral than this amount occur in the outer <1% of 
the radius (router = 1.454 x 10^^ cm for the 40kK benchmark 
model, as opposed to 1.463x10^^ cm for the average radial 
distance). In real H ll regions, dynamical effects dominate in 
these outer regions, where ionized gas is flowing away from 
the ionization front. With 65 cells on a side and the star at 
the centre, each cell has a width of 3% of the radius. Thus, 
the region we disregard (with hydrogen being J525% neutral) 
is much less than the width of each cell. 

We want to identify "extreme" predictions, defined by 
being outside of the range sh own by the other models. One 
way of describing them ( Pe auignot et alJl200 jl is the "isola- 
tion factor" , the ratio of our value to the next most extreme. 

For the 40kK benchmark, our code provides three ex- 
treme values, two within 4% of the nearest value. Our [S III] 
A9532+9069 prediction is 12.5% larger than that from the 
Rubin code (RR in the table). Rubin's code (RR), MOCASSIN 
(BE), and Harrington's (PH) are the only ones that use de- 
tailed radiative transfer rather than an approximation such 
as "outward only." 

For the 20kK benchmark, the results are very similar. 
We are extreme in six quantities, by <5% in all but two. 
The worst are [O III] A5949,5007 and (52/im + 88^im), (10% 
and 8% more than T. Kalman's xstar). We did not account 
for optical thickness in these lines, so our high predictions 
are understandable. Also, for this cool model our photon 
statistics for O"*"^ are problematic, as discussed above. Since 
ID codes can provide excellent spatial resolution by inserting 
zones wherever there are steep gradients, especially at the 
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outer edges where the temperature is changing rapidly, we 
feel that the benchmark results are satisfactory. 

One simple check is available for the cool model. There 
are relatively few He-ionizing stellar photons, and it is sim- 
ple to show that He absorptions strongly dominate H for He- 
ionizing photons. Each He ionizati on is followed by t he emis- 
sion of 0.96 H-ionizing photons l|Osterbroc3[l98i). Thus, 
the ratio of {He''')/(H^) is the ratio of the He- ionizing to 
H-ionizing stellar luminosities, 0.49 for a T = 20 kK black- 
body. Our code predicts 0.47, suggesting that our relatively 
coarse cell size is satisfactory for predicting He^/H^ and 
other ionic abundances. 

8.2 Ionization in a Shadow Region 

The structure of the shadow zone behind a small high- 
density clump of gas, such as those found in the plane- 
tary nebula NGC 7293, is of direct astrophysical interest. 
There have been several treatments of such a shadowed zone. 
Approximate static models have also be e n con sidered by 
IVan Blerkom fc Arnvl Jl972l) and 'Mathis (wid). An ana- 
lytical treatment is given bv lCanto et. al., (.1998i'l . who con- 
sidered dynamical evolution of the shadowed region and its 
ionization and temperature structure. These assume that 
the radiation field incident upon the shadow is the diffuse 
field within the H ll region, estimated by assuming that the 
nebula is thick to its own diffuse radiation and that there is 
no diffuse radiation produced by He. 

We tested the influence of He upon the ionization and 
temperature of a shadowed zone by comparing two models, 
one with the benchmark composition and one with the same 
H and heavy elements but no helium . Both used the geom- 
etry similar to that in ICanto et. alJ il998ft : a flat circular 
clump at 2; = 0, completely opaque to incident plane-parallel 
starlight out to its outer edge, so that it acts as a circular 
plate blocking stellar radiation incident normally upon it. 
We chose the uniform nebular density (100 cm~^) so that 
the stellar radiation is completely absorbed along the 6.2 pc 
z axis of the cube, even far from the clump. The ionizing 
spectrum is that of a 40000 K blackbody with an ionizing 
flux of 6.4 X 10^" cm~^ s~^. The radius of the occulting disk 
is 1.14 pc. We used periodic boundary conditions, meaning 
that when a photon exits one (x, y) face of the cube it reap- 
pears on the other side. Photons incident upon the z — 
face from within the cube are reflected back into the cube. 

The ionization in the shadow is indicated by Fig. |3] 
for the standard benchmark composition. We see that the 
diffuse radiation can penetrate relatively deeply into the 
shadow near the bottom of the figure (near the plate); the 
incident intensity is large because H is highly ionized. In 
this case, diffuse photons produced within relatively large 
distances can reach the shadow. At large z, near the top of 
the figure, the opacity in the directly illuminated regions is 
relatively large, the nebular radiation into the shadow is cor- 
respondingly low, and the ionized region has a sharp edge. 

Figure |1| shows T{r,z) vs. r/i?cioud, the distance from 
the axis of the shadow relative to the cloud's radius, and 
also T(r, z) vs. z, the distance along the shadow axis. Heavy 
lines are for standard composition; light lines show the no- 
He models. 

Fig. |lja) shows T{r,z), plotted against z, at r/_Rcioud 
= 1.17 (solid line), 1.01 (long dashes), and 0.964 (short 



dashes). The lines terminate when H becomes >95% neutral. 
We see that for the inner no-He model (light short dashes), 
T{r / Rcioud = 0.964) is much lower than the other tempera- 
tures at all z, as predicted in Cajito_ct^ al. (1998). H recom- 
bination radiation has a mean energy of only ~ 0.68fcsre — 
0.4 eV^. Photons released during the cascade following He 
recombinations (mostly the 19.8 eV 2^S — > l^S line) can re- 
lease ~ (20—13.6) eV = 6.3 eV in photoelectric heating with 
each ionization and penetrate ~3 times as far as H recombi- 
nations in the shadow. Thus, the He recombination photons 
have a major effect on the heating even though there are 
10 times as many from H. Even just outside of the shadow 
(the long dashed curves) the shadow temperature is about 
the same as far from the cloud (the solid curve), until the 
He becomes neutral in the gas fully illuminated by the star, 
at the point of the slight blip at z ~3.9 pc in the heavy 
solid curve. Outside of the shadow, we see the expected in- 
crease of T(r, z) with z, showing the effect of hardening the 
stellar radiation, for both helium- and no-helium models. 
(Since we did not account for cooling by coUisional excita- 
tion of H lines, the maximum temperature reached by the 
no-He model is probably significantly too large.) Within the 
shadow (short dashes), there is little change of T with z 
because the spectrum of the ionizing recombination radia- 
tion changes little along the axis of the shadow. Within the 
shadow, T{z) for the He models (heavy short dashes) is sig- 
nificantly less than outside, so the coUisionally excited line 
strengths are appreciably weaker. 

Fig.|ljb) shows T{r,z) plotted against (r/i?cioud) at two 
values of z: solid lines are for z = 0.1 pc; dashed for 2.6 pc. 
As before, the thick lines are for models with He, the thin for 
no helium. The inward drop of T{r) at r = i?cioud is obvious. 
We see that the low Te of the no-He models holds at all r. 
There is an increase in T(r) near the center for the helium 
models because the outer part of the shadow is mainly ion- 
ized by the soft recombination photons from H, while the He 
photons penetrate and deposit more energy when they are 
absorbed. There is no similar effect for the no-He shadow; all 
ionizing photons are soft. We see the penetrating power of 
the He photons into the shadow. Outside the shadow T{r) 
decreases somewhat with increasing r because the harder 
diffuse photons can reach the shadow more easily than the 
softer. 

This example shows the importance of diffuse helium 
radiation in shadow regions and may also help explain 
why He is less ionized t han H in the Diffuse Ionized Gas 
jRevnolds fc Tuftjll995^ . It shows that any treatment of 
the nebular diffuse radiation field inpinging onto shadowed 
regions must include a careful treatment of the helium re- 
combination radiation. We will investigate this effect and 
other 3D photoionization models in future papers. 

An anonymous referee provided helpful comments. JSM 
thanks a PPARC Visitors Grant to the University of St. An- 
drews, KW acknowledges support from a PPARC Advanced 
Fellowship. We thank Kirk Korista for comments and pro- 
viding the CLOUDY benchmark models and Mike Barlow and 



^ In common notation, the rate of kinetic energy loss per volume 
from recombinations is ne n(H"'") /JekTo, and the mean number 
of photons produced is rie n(H +) qr. Values of these functions 
are given m lOsterbrocd <1989l) . 
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Table 1. Lexington H ii Region Benchmark Results; H1I40" 



Line/H/3 


GF 


HN 


DP 


TK 


PH 


RS 


RR 


BE 


WME 




2.06 


2.02 


2.02 


2.10 


2.05 


2.07 


2.05 


2.02 


2.01 




0.119 


0.112 


0.113 


0.116 


0.118 


0.116 




0.114 


0.114 


P TT 9'i9'^-l- 


0.157 


0.141 


0.139 


0.110 


0.166 


0.096 


0.178 


0.148 


0.172 


P TTT 1 QOT-l-l QOQ 


0.071 


0.076 


0.069 


0.091 


0.060 


0.066 


0.074 


0.041 


0.078 


N TT 1 99//TY1 


0.027 


0.037 


0.034 




0.032 


0.035 


0.030 


0.036 


0.031 




0.669 


0.817 


0.725 


0.69 


0.736 


0.723 


0.807 


0.852 


0.742 




nn^n 

o.uooo 


0.0054 


n nn^n 

u. uuou 




U. UUU1 


u . uuou 


u . uuuo 


0.0061 


0.0057 


]VJ TTT "^7 '^l/m 
111 O 1 .OfJjLll 


0.306 


0.261 


0.311 




0.292 


0.273 


0.301 


0.223 


0.308 


W 1 oouut^uooo 


n nnQ4 
u.uuy^ 


u.uuoo 


u.uuoo 




u.uuoy 


n 0070 




oofi^ 


01 1 
U.Ul J. 




0.029 


0.030 


0.031 


0.023 


0.032 


0.024 


0.036 


0.025 


0.033 




1.94 


2.17 


2.12 


1.6 


2.19 


1.88 


2.26 


1.92 


2.23 


O TTT '^9-l-SK//m 


2.35 


2.10 


2.26 


2.17 


2.34 


2.29 


2.34 


2.28 


2.42 


O TTT '^nn7-l-ZLQ'^Q 


2.21 


2.38 


2.20 


3.27 


1.93 


2.17 


2.08 


1.64 


2.34 


O in 4363 


0.00235 


0.0046 


0.0041 


0.0070 


0.0032 


0.0040 


0.0035 


0.0022 


0.0044 


Nc II 12.8^01 


0.177 


0.195 


0.192 




0.181 


0.217 


0.196 


0.212 


0.192 


Nc III 15.5/^m 


0.294 


0.264 


0.270 


0.35 


0.429 


0.350 


0.417 


0.267 


0.263 


Ne III 3869+3968 


0.084 


0.087 


0.071 


0.092 


0.087 


0.083 


0.086 


0.053 


0.063 


S II 6716+6731 


0.137 


0.166 


0.153 


0.315 


0.155 


0.133 


0.130 


0.141 


0.18 


S II 4068+4076 


0.0093 


0.0090 


0.0100 


0.026 


0.0070 


0.005 


0.0060 


0.0060 


0.0094 


S III 18.7/im 


0.627 


0.750 


0.726 


0.535 


0.556 


0.567 


0.580 


0.574 


0.623 


S III 9532+9069 


1.13 


1.19 


1.16 


1.25 


1.23 


1.25 


1.28 


1.21 


1.44 


10^ X A(BC3645)/A 


4.88 




4.95 




5.00 


5.70 




5.47 


4.96 


^inncr 


7719 


7668 


7663 


8318 


7440 


7644 


7399 


7370 


7743 


< T[NpNe] >/K 


7940 


7936 


8082 


8199 


8030 


8022 


8060 


7720 


8183 


Rout/lOl^cm 


1.46 


1.46 


1.46 


1.45 


1.46 


1.47 


1.46 


1.46 


1.45 


< Hc+H+ > 


0.787 


0.727 


0.754 


0.77 


0.764 


0.804 


0.829 


0.715 


0.771 



°GF, Fcrlaiid's CLOUDY; HN, H. Nctzer's ION; DP, D. Pequinot's NEBU; TK, T. Kallman's XSTAR; 
PH, J. P. Harrington's code; RS, R. Sutherland's mappings; RR, R. Rubin's nebula; 
BE, B. Ercolano's mocassin; WME, this paper. 



Table 2. Lexington H ii Region Benchmark Results: HII20 



Line/H/3 


GF 


HN 


DP 


TK 


PH 


RS 


RR 


BE 


WME 


H^/1036ergs-i 


4.85 


4.85 


4.83 


4.90 


4.93 


5.04 


4.89 


4.97 


4.87 


He I 5876 


0.0072 


0.008 


0.0073 


0.008 


0.0074 


0.0110 




0.0069 


0.00684 


C II 2325+ 


0.054 


0.047 


0.046 


0.040 


0.060 


0.038 


0.063 


0.042 


0.053 


N II 122^jrn 


0.068 




0.072 


0.007 


0.072 


0.071 


0.071 


0.071 


0.067 


N II 6584+6548 


0.745 


0.786 


0.785 


0.925 


0.843 


0.803 


0.915 


0.846 


0.778 


N II 5755 


0.0028 


0.0024 


0.0023 


0.0029 


0.0033 


0.0030 


0.0033 


0.0025 


0.0025 


N III 57.3/irn 


0.0040 


0.0030 


0.0032 


0.0047 


0.0031 


0.0020 


0.0022 


0.0019 


0.0049 


O I 6300+6363 


0.0080 


0.0060 


0.0063 


0.0059 


0.0047 


0.0050 




0.0088 


0.055 


O II 7320+7330 


0.0087 


0.0085 


0.0089 


0.0037 


0.0103 


0.0080 


0.0100 


0.0064 


0.0089 


O II 3726+3729 


1.01 


1.13 


1.10 


1.04 


1.22 


1.08 


1.17 


0.909 


1.12 


O III 52+88/im 


0.0030 


0.0026 


0.0026 


0.0040 


0.0037 


0.0020 


0.0017 


0.0022 


0.0044 


O III 5007+4959 


0.0021 


0.0016 


0.0015 


0.0024 


0.0014 


0.0010 


0.0010 


0.0011 


0.0026 


Nc II 12.8Mm 


0.264 


0.267 


0.276 


0.27 


0.271 


0.286 


0.290 


0.295 


0.289 


S II 6716+6731 
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Figure 1. HII40 benchmark, showing our Monte Ccirlo results (crosses), mocassin (squares), and the cloudy model (solid line). 



Barbara Whitney for comments on early versions of the pa- 
per. The MOCASSIN tests were carried out on the Enigma 
SunFirc Cluster, at the HiPorSPACE Computing Centre, 
UCL, which is funded by PPARC. 
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Figure 2. HII20 benchmark, showing our Monte Carlo results (crosses), MOCASSIN (squares), and the CLOUDY model (solid line). 



1996, ApJ, 462, 326 
Hummer, D. G., & Storey, P.J. 1998, MNRAS, 297, 1073 
Hummer, D. G. 1994, MNRAS, 268, 109 
Katz, N., Weinberg, D. H., & Hernquist, L. 1996, ApJS, 

105, 19 

Kingdon, J. B., & Ferland, G. J. 1996, ApJS, 106, 205 
Lauroesch, J. T., Meyer, D. M., & Blades, J. C. 2000, ApJ, 

543, L43 
Lucy, L. B., 2003, A&A, 403, 261 
Lucy, L. B., 1999, A&A, 344, 282 
Mathis, J. S. 2000, ApJ, 544, 347 
Mathis, J. S., 1976,ApJ, 207, 442 

Mathis, J. S., Whitney, B. & Wood, K. 2002, ApJ, 574, 812 
Miller, W. W., HI, & Cox, D. R1993, ApJ, 417, 579 
Minter, A. H., & Spangler, S. R., 1996, ApJ, 458, 194 
Nussbaumer, H., & Storey, P. J. 1983, A&A, 126, 75 
Nussbaumer, H., & Storey, P. J. 1984, A&AS, 56, 293 
Nussbaumer, H., & Storey, P. J. 1987, A&AS, 69, 123 
Och, S. R., Lucy, L. B., & Rosa, M. R. 1998, A&A, 336, 
301 

Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebu- 
lae and Active Galactic Nuclei, University Science Books, 
Mill Valley, CA 

Pequignot, D., et al. 2001, in "Spectroscopic Challenges of 
Photoionized Plasmas", G. Ferland & D. W. Savin, eds.. 



ASP Conference Series Vol. 247. San Francisco: Astro- 
nomical Society of the Pacific, 533 

Pradhan, A. K., & Peng, J. 1995, in Space Telescope Sci- 
ence Institute Symposium Series No. 8, Ed: R.E.Williams 
and M.Livio, Cambridge University Press, p8 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flan- 
nery, B. P. 1992, Numerical Recipes in C-l-l- (2d ed.; New 
York: Cambridge Univ. Press) 

Reynolds, R. J., & Tufte, S. L. 1995, ApJ, 448, 715 

Reynolds, R. J., Haffner, L. M., & Madsen, G. J., 2002, in 
"Galaxies: the Third Dimension", M. Rosada, L. Binette, 
& L. Arias, eds. (Astr. Soc. Pacific,: San Francisco) 

Reynolds, R. J., Haffner, L. M., & Tufte, S.L. 1999, ApJ, 
525, L21 

Savage, B. D., & Sembach, K. R. 1996, ARAA, 34, 279 
Seaton, M. J. 1959, MNRAS, 119, 81 

Simon, R., Jackson, J. M., Clemens, D. P., & Bania, T. M. 

2001, ApJ, 551, 747 
Spaans, M., 1996, A&A, 307, 271 

Spitzer, L. 1978, Physical Processes in the Interstellar 

Medium, New York, Wiley 
Stanimirovic, S. & Lazarian, A. 2001, ApJ, 551, L53 
Stasiriska, G., 2002, lectures at the XIII Canary Islands 

Winterschool on Cosmochemistry, astro-ph/0207500 
Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41 



12 Wood, Mathis, & Ercolano 




Figure 3. The ionization structure, H'^/H, of a shadowed region produced by an impervious circular cloud, with ionizing stellar radiation 
incident from below as described in §8.2. 
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Figure 4. (a) The nebular temperature, T{r,z), within a region shadowed by a circular cloud that acts as an opaque plate, as shown 
in Fig. 3. Stellar photons are injected parallel to the z axis (vertical in Fig. 3); r is the radial distance from the axis of the shadow. 
The radius of the cloud is 1.14 pc. Heavy lines show the case of the benchmark nebular composition, roughly solar; light lines are solar 
except that there is no helium. Solid lines; r/i?(cloud) = 1.17; long dashes:, r/R(cloud) = 1.01; short dashes:, r/-R(cloud) = 0.96. The 
figure shows the dramatic effects that He recombination radiation has on the physical conditions within the shadow (see text), most 
importantly, the very low Te if there is no He. (b) T{r,z) as a function of z at r/fl(cloud) = 0.095 pc (solid lines) and z = 2.6 pc 
(dashed). 



